The static equilibrium design is one of the most fundamental studies in tokamak fusion science. Experimental Advanced Superconductive Tokamak (EAST) is a D-shaped cross section tokamak with 14 independent full superconductive poloidal field (PF) coils (PF7 and PF9, PF8 and PF10 are connected in series, respectively) which are symmetrically located outside the vessel, as shown in Fig. 1. The plasma current, shape, and position are modified by electromagnetic fields generated by PF coils.^[1] In order to reduce the heat load on the divertor plates, the snowflake (SF) divertor configuration was proposed and studied.^[2] By now, the SF has been established in TCV, NSTX, and DIII-D tokamak,^[3] but in EAST tokamak, only quasi-snowflake (QSF) configuration could be achieved because of the limitation of PF coils, and QSF discharge has been carried out in EAST tokamak since 2014.^[4]

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. Cross section of EAST and plasma.

The simulation and experimental results of QSF discharges show that the poloidal flux distribution is not steady around the X-point region.^[5] Additionally, the formation of QSF configuration with lower X-point needs high currents of PF6 and PF14, which may exceed the limitation of the coils. QSF configuration design and control remain challenging. Various codes and numerical models have been developed to optimize the EAST plasma equilibrium.^[6,7] EFIT code is one of the most widely used plasma boundary reconstruction codes, which is described in Ref. [8]. Tokamak simulation code (TSC), which is used to analyze experiment data off-line, performs the numerical simulations of the nonlinear evolution of plasma and the magnetic field. In Ref. [9], the EAST digital model is constructed and QSF discharge is reproduced by TSC using RZIP method. CREATE-NL code solves the fixed-boundary dynamic plasma equilibrium problem based on finite element method.^[10] It is indicated that the design of QSF configuration consists of two steps: firstly creating a configuration with a second null point based on a standard single null (SN) divertor configuration, which have almost the same plasma boundary, and secondly modifying the plasma shape and PF current distribution while satisfying the technical constraints.

In the EAST QSF configuration design, it is found that the desired PF6 and PF14 currents usually exceed the PF current limitation (12 kA/turn).^[5] Obviously, this is not suitable. By modifying the plasma shape, the PF currents can be changed. Thus, in this paper, optimization of plasma shape of static QSF equilibrium using the orthogonal method is proposed with the intention of reducing the currents of PF6 and PF14 as small as possible. A fixed boundary equilibrium solver based on a non-rigid plasma model is used to calculate the PF current distribution at first. Then the QSF configuration with total plasma current 300 kA is designed. According to the orthogonal method, the currents of PF6 and PF14 are regarded as the target functions, while the shape parameters triangularity, elongation, and squareness are regarded as variables. Equilibrium calculations under different variable combinations are carried out. Based on the results, the plasma shape is redefined. The simulation shows that the target PF currents of optimized QSF equilibrium are reduced.

Shoot 56573 is a typical QSF discharge with 250 kA and toroidal field 1.8 T. The reconstructed plasma equilibrium at 5 s by off-line EFIT code is set as the reference equilibrium. The plasma configuration is shown in Fig. 2, and the equilibrium parameters are shown in Table 1. Based on the reference equilibrium, the QSF configuration with 300 kA is designed, which will be discussed in the next section.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. Reconstructed plasma equilibrium of shoot 56573 at 5 s by off-line EFIT code. The blue solid line denotes the plasma boundary and the gray dashed lines denote poloidal flux surfaces. The top-outer squareness is defined by sqto = CD/BD, where D is the middle point of AB.

Table 1.

Table 1.

Table 1. Equilibrium parameters of Shoot 56573 at 5 s. . Parameters Symbol Specification Plasma current/kA Ip 250 Toroidal field/T Bt 1.8 Poloidal beta βp 0.6 Internal inductance li 1.6 Major radius/m R0 1.9 Minor radius/m a0 0.45 Elongation e 1.34 Electron density/m−3 ne 2.25 × 1019

Table 1. Equilibrium parameters of Shoot 56573 at 5 s. .

In the discharge of EAST tokamak, the PF coils are used for both ohmic heating and plasma shaping. Each PF coil is powered by an individual power supply and the desired current in the coil is composed of a feedforward trajectory value and a feedback-calculated value.^[11] The feedforward value for a new plasma configuration, which is usually proposed based on one produced in the experiment, is calculated by solving fixed-boundary equilibrium problem.

3.1. Fixed boundary equilibrium solver

In the axisymmetric tokamak configuration, plasma equilibrium is described by the Grad–Shafranov (GS) equation derived from the steady state ideal magneto-hydrodynamic (MHD) equations,^[12] and is given by

where j_φ is toroidal component of the plasma current density distribution in cylindrical coordinates (R, Z, Φ), μ₀ is the magnetic permeability of free space, ψ is the normalized magnetic flux, P(ψ) is plasma pressure, and F(ψ) is the poloidal current function defined by F(ψ) = rB_ϕ where B_ϕ is the toroidal component of the magnetic field.

According to the GAQ plasma model proposed by Brown,^[13] j_φ can be described as

Here, , where ψ_max is the value of ψ on the magnetic axis and ψ_lim is the value on the plasma boundary. The parameter λ is a normalization factor specified in terms of the total current I_p. The quantity r₀ is major radius of the plasma. The quantities α_m, α_n, γ_m, and γ_n need to be adjusted in order to match the internal inductance l_i, which describes the peakedness of the current profile. The parameter β₀ is a scalar quantity related to the poloidal beta β_p.

For fixed plasma boundary and current profile, the calculation of PF current distribution essentially is a process to minimize the cost function f(I_c, ψ),^[14] which is defined as

where I_c = (I₁,…, I_n) is the vector of PF coil currents, ψ(M) is the flux at the control points of plasma boundary, Γ_d is the desired plasma boundary, Ω_X is the region near the X-point, w_M, w_X, and w_i are weighting coefficients specified by user.^[15]

It is derived that the equation (3) can be converted into an overdetermined equation system, which is composed of a group of least square solutions, and it is given by

where M_lc is mutual inductance between PF coil currents and control points, G_xrc and G_xzc are Green’s functions between magnetic field (horizontal and vertical) with PF coil currents, I is the identity matrix, ψ₀ is the constraint value of flux on the plasma boundary, is the flux at the control points produced by plasma current, and are horizontal and vertical magnetic field at the X-points produced by plasma current.

In order to solve the above equations, a fixed boundary equilibrium solver which contains the current profile adjustment module has been developed in MATLAB software (MathWorks, USA) based on the GAQ plasma model. The flow chart of the computation is shown in Fig. 3. To verify the developed solver, the equilibrium of shoot 56573 at 5 s is calculated based on the EFIT reconstructed plasma shape. The plasma model is set up using Eq. (2), and the values of α_m, α_n, γ_m, γ_n, and β₀ are adjusted to obtain the same l_i and β_p with EFIT plasma model. The equilibrium calculated by the proposed solver is shown in Fig. 4, and the PF current distribution is compared with that of EFIT code. The result indicates that the error of calculated PF currents between two codes is less than 3%, which proves the reliability of the proposed solver.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. Flow chart of the orthogonal method.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. Plasma equilibrium of shoot 56573 at 5 s calculated by the proposed solver. PF currents calculated by the proposed solver and EFIT code are compared.

By exploiting the solver, the plasma equilibrium mentioned in Section 2 is modified. The target plasma configuration has the same l_i, β_p, and plasma shape with the reference equilibrium, except that the I_p is set to 300 kA. The result of equilibrium calculation is shown in Fig. 5, which reveals that the currents of PF6 and PF14 exceed the maximum limited value 12 kA. To solve this problem, the equilibrium parameters should be modified. However, the full optimization of equilibrium is very complicated. This paper focuses on the plasma shape modification, more specifically, tweaking the shape parameters elongation e, bottom triangularity dbot, and squareness. The squareness is divided into top-inner squareness sqti, top-outer squareness sqto, bottom-inner squareness sqbi, and bottom-outer squareness sqbo. The optimized equilibrium with the redefined plasma shape is supposed to have the PF current distribution within limits.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. Plasma configuration before optimization. The currents of PF6 and PF14 exceed the maximum limited value 12 kA.

3.2. Orthogonal method

The orthogonal method,^[16] based on the combinatorial theory, is an efficient way to test pair-wise interactions. It reduces the number of tests by choosing representative variable combinations that covers the whole test areas. With this method, the parameters that affect the test results are defined as “factors”, and the different values of factors are “levels”. Depending on the number of the factors and levels, the orthogonal array L_m(p^q) can be constructed. Here, m is the total number of simulations, p and q are the numbers of levels and factors, respectively.

To optimize the plasma shape, the currents of PF6 and PF14 are set as two independent target functions, and the shape parameters elongation, triangularity, and squareness are considered as factors. The levels are decided with the following constrains. (i) The number of the levels are within an appropriate range. Too many levels will increase the computational complexity, while inadequate level may misread the change trend of the factor. (ii) The plasma boundary has a minimum of 40 mm clearance gap from the first wall. (iii) The strike points are on the vertical targets. Given the above, five levels for each factor are selected and a orthogonal table as demonstrated in Table 2 is constructed. For the analysis of the results of the orthogonal tests, two important indexes need to be calculated: the average value k_i and the range R. The average value is used to assess the influence of specific factor in its level on the target functions, which is defined as

where i is the level index, J_in is the simulation result for specific factor at level i, and N_i is the number of simulations of the corresponding factor at level i.

The range R is used to evaluate the influence weight of different factors (the larger the range R is, the higher the influence weight is), which is defined as the difference between maximum and minimum k_i for the corresponding factor:

Table 2.

Table 2.

Table 2. Orthogonal array with 6 factors and 5 levels. . Column 1 2 3 4 5 6 PF6 PF14 Factor e dbot sqto sqti sqbo sqbi current/kA current/kA 1 1.32 0.29 0.39 0.44 0.25 0.21 8.86 14.9 2 1.32 0.30 0.40 0.45 0.26 0.22 9.19 14.2 3 1.32 0.31 0.41 0.46 0.27 0.23 9.9 13.8 4 1.32 0.32 0.42 0.47 0.28 0.24 10.6 12.9 5 1.32 0.33 0.43 0.48 0.29 0.25 11.4 11.7 6 1.33 0.29 0.40 0.46 0.28 0.25 9.9 12.4 7 1.33 0.30 0.41 0.47 0.29 0.21 8.4 11.8 8 1.33 0.31 0.42 0.48 0.25 0.22 10.2 11.4 9 1.33 0.32 0.43 0.44 0.26 0.23 9.6 11.4 10 1.33 0.33 0.39 0.45 0.27 0.24 11.4 12.5 · · · · · · · · · · · · · · · · · · · · · · · · · · · 25 1.37 0.33 0.42 0.46 0.26 0.21 8.1 9.3 k1 13.82 13.01 13.50 12.53 11.16 12.6 k2 12.02 12.58 13.03 12.6 11.18 12.42 k3 11.87 11.38 11.60 12.9 11.30 12.36 k4 11.43 11.22 11.35 12.29 12.07 12.21 k5 11.10 11.03 11.08 12.2 12.20 12.18 R 2.72 1.98 2.42 0.7 1.04 0.42

Table 2. Orthogonal array with 6 factors and 5 levels. .

According to Table 2, 25 simulations are conducted and the results are shown in Table 2. According to Eqs. (5) and (6), the average value of each level and the range of each factor are calculated, as shown in Fig. 6. Table 2 only shows the case that target function is PF14 current for brevity. It is clear that the current magnitude of PF6 increases monotonically with the increment of sqbo and sqbi, and for factors e, dbot, sqto, and sqti, the PF6 current reaches the minimum when they are 1.36, 0.33, 0.46, and 0.48, respectively. The current magnitude of PF14 decreases with the increment of e, dbot, sqto, and sqbi. For sqti and sqbo, the PF14 current reaches the minimum when they are 0.48 and 0.25, respectively. Furthermore, the result indicates that among all these factors, sqbo and e have the biggest impact on PF6 current and PF14 current.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. The dependence of PF6 and PF14 currents on plasma shape parameters.

With the currents of PF6 and PF14 considered together, the best shape parameters combination is optimized, as shown in Table 3. Finally, the equilibrium with redefined plasma shape is computed, as shown in Fig. 7. It can be seen that the currents of PF6 and PF14 are reduced to 8.417 kA and −10.956 kA, respectively.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. Plasma configuration after optimization. The currents of all PF coils are within allowable range.

Table 3.

Table 3.

Table 3. Parameters comparison. . Variables Preliminary design Final design (Before optimization) (After optimization) R0/m 1.9 1.9 a0/m 0.45 0.45 e 1.34 1.36 dbot 0.31 0.33 sqto 0.41 0.43 sqti 0.46 0.48 sqbo 0.27 0.25 sqbi 0.21 0.20 PF6/kA 12.009 8.417 PF14/kA –13.729 –10.956

Table 3. Parameters comparison. .

The main objective of the QSF configuration is to increase the flux expansion and thus reduce the heat flux density to the target. From this point of view, the flux expansion of the QSF configuration before and after optimization are compared. Figure 8 shows the scrape-off layer (SOL) geometry near the null point, and the flux surfaces near the separatrix in both cases are assumed to be situated at the distance of 1.75 cm from the separatrix at the equatorial plane. The flux expansion of the optimized configuration is 7.15 cm, which is larger than 5.03 cm before optimization. Thus the heat flux density to the target can be reduced after optimization.

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. The structure of magnetic flux surfaces near the null-point of the QSF configuration.

By exploiting the proposed optimization method, we also designed the QSF equilibrium with upper X-points and I_p = 300 kA. Unlike the QSF equilibrium with lower X-points, the QSF configuration with upper X-points needs high PF5 and PF13 currents, which may exceed the PF current limitation.^[5] Thus, based on a typical QSF equilibrium (shoot 71543 at 3 s) with I_p = 250 kA, the new equilibrium is designed with PF5 and PF13 currents set as target functions and elongation, top triangularity dtop, and squareness set as variables. The design result is shown in Table 4. Before optimization, the desired PF5 and PF13 currents exceed the PF current limitation, and all PF currents after optimization are within allowable range.

Table 4.

Table 4.

Table 4. Design result of QSF configuration. . Variables Reference Preliminary design Final design Configuration (Before optimization) (After optimization) Ip/kA 250 300 300 R0/m 1.85 1.85 1.85 a0/m 0.42 0.42 0.42 e 1.95 1.95 2.0 dtop 0.54 0.54 0.56 sqto 0.35 0.35 0.32 sqti 0.19 0.19 0.2 sqbo 0.37 0.37 0.39 sqbi 0.41 0.41 0.43 PF1/kA 1.230 2.013 –1.798 PF2/kA 6.745 7.354 5.904 PF3/kA 6.392 7.035 6.153 PF4/kA 5.189 6.723 5.448 PF5/kA 9.657 12.405 9.684 PF6/kA 6.335 5.354 4.035 PF7/9/kA 3.514 4.653 2.955 PF8/10/kA 4.089 5.094 2.856 PF11/kA 3.527 3.267 2.921 PF12/kA 0.298 1.675 2.915 PF13/kA –9.235 –14.256 –10.956 PF14/kA –1.763 –4.567 –5.286

Table 4. Design result of QSF configuration. .

In order to validate the designed equilibrium, an experiment was conducted as shoot 73559 in EAST. The preprogrammed PF currents at 3 s of shoot 73559 were set to the designed value in Table 4. Figure 9 shows the measured experimental data of PF currents and total plasma current, which are basically consistent with the designed values at 3 s. Figure 10 gives the EFIT reconstructed plasma equilibrium at 3 s, and the plasma shape is compared to the designed shape, which shows a good consistency.

Fig. 9.

	Figure Option ViewDownloadNew Window
	Fig. 9. Experimental data of PF currents and total plasma current of shoot 73559.

Fig. 10.

	Figure Option ViewDownloadNew Window
	Fig. 10. EFIT reconstructed plasma equilibrium of shoot 73559 at 3 s. The designed plasma boundary (red line) is superimposed on the EFIT boundary. The black dots represent the reference plasma shape.

This paper studies the static QSF equilibrium of EAST tokamak. Generally, the formation of the second X-point needs large PF5/6 and PF13/14 currents. Thus it is important to take into account the restraint of coil currents during equilibrium design. A fixed boundary equilibrium solver has been developed to calculate the plasma equilibrium, and the equilibrium of shoot 56573 at 5 s is computed by the developed solver and EFIT code, respectively. The plasma models of two codes have different analytical forms. Therefore, the plasma current profiles calculated by the two codes are not completely the same, and the error of the calculated PF currents between two codes is less than 3% according to the calculation results. Then the influence of plasma shape on PF coil current distribution is analyzed by orthogonal method as well as the optimization of plasma shape. The range analysis indicates that the parameters impacting the PF6 current from high to low are sqbo, dbot, sqbi, sqti, sqto, e, while for PF14 current, the rank is e, sqto, dbot, sqbo, sqti, sqbi. With shape parameters e, dtop, sqto, sqti, sqbo, and sqbi changed from 1.34, 0.31, 0.41, 0.46, 0.27, and 0.21 to 1.36, 0.33, 0.43, 0.48, 0.25, and 0.20, respectively, the current magnitude of PF6 and PF14 are reduced from 12.009 kA and 13.729 kA to 8.417 kA and 10.956 kA, respectively. Finally, the experiment is conducted to validate the optimization of the QSF with upper X-points, and the desired QSF plasma configuration is obtained with the preprogrammed PF currents.